Proportion Regulation in Globally Coupled Nonlinear Systems 



Tsuyoshi Mizuguchi * and Masaki Sano ' 
Research Institute of Electrical Communication 
Tohoku University, Sendai 980-11, Japan 

Abstract 

As a model of proportion regulation in differentiation process of biological 

system, globally coupled activator-inhibitor systems are studied. Formation 

and destabilization of one and two cluster state are predicted analytically. 

Numerical simulations show that the proportion of units of clusters is chosen 

within a finite range and it is selected depend on the initial condition. 
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The regulation of proportion among different cell types in a tissue is a general and 
important aspect of biological development. It is well known that the proportion between 
the two different cell types is roughly constant irrespective of the slug size of cellular slime 
mold Dictyostelium discoideum (Dd) amoebae |I|-fJJ. Initially the same type of aggregative 
cells, when dissociated, randomly mixed, and reaggregated, differentiate into two types 
cells (prespore and prestalk cells) without pattern formation. It is known now that cell 
differentiation starts independently on the cell position, and later cell sorting forms the two- 
zoned prestalk-prespore pattern in slug of Dd Similar regulation mechanism can be 
observed in caste populations of social insects such as ants and bees In the division 
of work, proportion is regulated irrespective of the size of society nor the artificial partial 
extinction by the experimenter. 

No theoretical model exists to describe the proportion regulation. Any pattern formation 
model such as Turing type instability with diffusive coupling || is incompatible with the 
observation that the Dd cells start to differentiate independent to their positions. A large 
population of identical units interacting equally to the other units (globally coupled nonlinear 
system) is a good candidate to describe the phenomena. It is an idealized model of the cases, 
when the diffusion length of chemical factor, e.g. differentiation inducing factor (DIF) or 
pheromone, is large enough compared to the cell size, or when the individual units moves 
around to interact with others. 

Recently, globally coupled chaotic map flJTCfl and globally coupled oscillators [PH T6 



are studied and interesting phenomena including clustering and their destabilization are 
observed. However, analysis of cluster state is difficult for these systems because the unit 
itself is complex enough. It is also doubtful that chaos or oscillation is playing essential roles 
in proportion regulation of biological system such as Dd. In this respect, a minimum model 
of clustering is preferable. 

Our model of globally coupled system is composed of N activator-inhibitor type units 
which have two variables u and v. The dynamics of each unit is modeled as 
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iij = au,j — bvj — u'j +Ki(u — Uj), 
Vj = cuj — dvj +K 2 (v — Vj), 

Y N j N 

u = — y^Mj and v = — Y^ik. 

1=1 i=i 

Here, two component Uj and Vj of j-th unit are considered as activator and inhibitor by 
assuming that a,b,c, and d are positive. Each unit couples with all other units through the 
averaged field u and v. We assumed that both K\ and K 2 are non-negative. They can be 
regarded as the susceptibility of each component because this type of global coupling can 
be considered as the fast limit of diffusion velocity. 

First, we investigate the properties of individual unit by setting K\ = K 2 = 0. Steady 
stationary solutions (uq, vq) are easily solved by setting ii — v — 0. Depend on the parameter 
s = ad — be, the number and the stability of the fixed point changed. Linear stability of 
these fixed points can be analyzed by setting u = uo + 5u,v = vq + 5v,\Su\,\Sv\ 1, and 
5u = 5u e xt ,5v = 5v e xt . Linearization of (0) leads the eigenvalue equation 

= A 2 - (a - d - 3w 2 )A - s + 3du 2 . (2) 

The fixed point (uo,v ) is stable if the conditions > a — d — 3mq and < — s + 3du^ are 
satisfies. We assume that a < d and s < from now. In this case, the trivial solution (0, 0) 



is the unique attractor of the dynamical system of one individual unit. |T7 . 

Let us consider a one cluster state which is defined as a state that every unit has the 
same value ("U(i), ^(i))) i>e. 

(uj,Vj) = (w(i),V(i)) for j=l,...,N. 

It is generally difficult to analyze the stability of cluster state in the globally coupled system 
because we must solve eigenvalue problem of 2N dimensional matrix. Following analyzing 
method, however, we obtain sufficient condition for destabilization of a cluster state easily. 
First, N units are assumed to form M clusters state {(«(»), v (*))}, i = 1,...,M. Next, we 
consider one additive unit (test unit) (u(t),v(t)) in that M clusters state and make an 



approximation that the effect to the test unit from N units is simply external force. This 
approximation is justified in the limit of iV — > oo. By investigating the stability of this test 
unit, we argued the stability condition of original M clusters state as follows. Noting that the 
test unit has at least M "entrained" solutions to each cluster, i.e., (u(t),v(t)) = f(i)), 
the linear stability of these entrained solutions can be analyzed. If one of the entrained 
solutions is unstable, we conclude that the original cluster state is unstable. We named this 
stability analysis method of the cluster state Test Unit Analysis (TUA) ||18|| . 

Now, we carry out TUA for the one cluster state {(0, 0)} given in the previous paragraph, 
i.e., we consider a stability of test unit (u(t),v(t)) in the external force created by the N 
units in one cluster state. In this case both of the average fields u and v vanish and equations 
for the test unit are written in the form: 

u = (a — Kt)u — bv — u 3 , 

(3) 

v = cu — (d + K 2 )v 

Fixed point of test unit can be obtained by setting ii — v — 0. Using v = cu/(d + K 2 ) given 
from (Q), u satisfies 

= h x [u) = {d + K 2 )u 3 + {(d + K 2 )(K 1 - a) + 6c) w. 

Note that hi(0) = because the test unit has the entrained solution (u,v) = (0,0). To 
investigate the stability of one cluster state, we analyze the test unit linear stability around 
entrained solution (0, 0). By setting u = 5u and v = Sv, linear stability analysis of (||) leads 
the eigenvalue equation: 

= A 2 - (a - d - K x - K 2 )\ - (a - K t ){d + K 2 ) + be. (4) 

The stability conditions of the entrained solution of the test unit to the one cluster state are 
now given as: 

0>a-d-K 1 -K 2 (5) 
0<-s-aK 2 +K 1 (d+K 2 ). (6) 
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From the condition for existence and stability of one cluster solution in the non-coupling 
case, (|5|) is automatically satisfied . Therefore the critical condition for the stability is given 
by R.H.S of (|6|) equals to where a pitch-fork bifurcation occurs. Although this stability 
condition is for the entrained solution of the test unit, it is obvious that the original one 
cluster solution is unstable if (^) is broken. From the fact mentioned above, one cluster 
state is linearly unstable when K 2 > K 2c . Fig. 1 shows a result of numerical simulation for 
K 2 < K 2c . Parameters are iV = 100, a = 0.4, b = 1, c = 0.5, d = 1, K x = 0, and K 2 = 0.2. 
Simple Euler method with dt = 0.01 are adopted. One cluster state (Fig. 1(c)) is realized 
from a uniform random initial condition (Fig. 1(b)). When K 2 > K 2c , the one cluster state 
becomes unstable and each unit separates into two subpopulations, i.e., 

. , ( w (i)> u (i)) for * = !> •••>%) 
{u h Vi) = < 

[ (u( 2 ),v {2) ) for i = N {1) + 1, N 

Here, Nm is a number of units which belongs to the first cluster. This state is defined as 
two clusters state and we focus on it in the next paragraph. 

For simplicity we assumed that K\ = 0. Let a two clusters state (w(i), f(i)),(u(2), U(2)) 
having its proportion p : 1 — p, here p is defined as p = N^/N and < p < 1 are satisfied. 
Because the averaged fields u and v become pum + (1 — p)ur 2 ) and pvm + (1 — p)vr 2 \, 
respectively, «(i), W(i), W(2), ^(2) satisfy the following equations: 

= au(i) — — 

= cu ( i) - dv (1) + K 2 (l - p)(v( 2 ) - V(i)) 

= aU( 2 ) - bV( 2 ) - W( 2 ) 

= CU( 2 ) - C?t)(2) + K 2 p(v {1) - tl( 2 )) 

Eliminating t>(x) and v (2), transforming from (ii(i), W(2)) to («(i), 0) where is a new variable 
defined by = u^/u^i), one can get two equations which are easily analyzed: 

2 _ s + aK 2 (l -p) - aK 2 (l —p)<f> 
U{1) ~ d + K 2 (l-p)-K 2 (l-p)<jfi 
(s + aK 2 p)(j) - aK 2 p 
~ (d + K 2 p)<p* - K 2 p ' 
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obeys 

(0 - l)(bcK 2 (l - p)<f - s(d + K 2 )<p 2 - s(d + K 2 )<p + bcK 2 p) = 0. 

Note that = 1 expresses Um = Uc 2 \, i.e., one cluster state. The solution must satisfies the 
inequality ufa > 0. Therefore the condition for existence of two clusters solution is : 

K 2 > K 2c = -s/a. (7) 

Under the condition (0), we use TUA again to analyze the stability of the two clusters 
state, i.e., we consider a stability of test unit entrained solution in the external force created 
by the N units in the two clusters state. The test unit equations are: 

u = au — bv — u 3 

v = cu- dv + K 2 (pv(i) + (1 - p)v(2) - v). 

By setting u = v = 0, we obtain fixed points of the test unit. Note that there exist two 
entrained solutions to the cluster (um,^)) and (up), v (2))- To investigate the stability of two 
cluster state, we analyze the test unit linear stability around entrained solution (um,vn)). 
By setting u = U(\) + 6u,v = V(\) + Sv, linearization around (w(i),f(i)) leads the eigenvalue 
equation: 

= A 2 - (a - d - K 2 - 3u 2 {1) )\ - (a - 3u 2 {1) ){d + K 2 ) + be. (8) 

The entrained solution becomes unstable if the constant term in (0) becomes 0. The stability 
condition is 

> {(aK 2 + s)((2ad - 3bc)K 2 + 2dsf(d + K 2 p)} 

x{-9bcK 2 p + 2{ad + 3bc)K 2 + 2ds}. (9) 

From the fact that first braces of R.H.S. of (^) are positive definite, the last braces determine 
the stability. Bifurcation line is given by solving about p, i.e., 

p = p c (K 2 ) = -L(^£ + 2 (ad + 36c)), (10) 
9oc K 2 
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where a transcritical bifurcation occurs. For example, p c = 7.6/9 — 0.2/4.5.^2 for the case 
a = 0.4,6 = l,c = 0.5, d = 1. Typical phase diagram is shown in Fig. 2. In the case 
of K% < K2C1 there does not exist two clusters solution. In the case of K2 > K2 C , on 
the other hand, one of the two clusters solution with a proportion lain in the region B is 
realized. The two clusters state in the region C is linearly unstable and is not realized. 
Note that the bifurcation diagram is symmetric to p = 0.5 because the proportion of the 
other cluster is 1 — p. Therefore possible proportion has a minimum p m j n = 1 — p c {K.2) 
and maximum p max = Pc{K-2) value for given K2. To investigate the dynamical process 
we perform numerical simulations. Fig. 3 shows the formation of two clusters state from a 
uniform random initial condition with K 2 = 0.3 > K 2c . Other parameters are the same as 
in Fig. 1. At T = 500, two clusters state is selected with a proportion p : 1 — p = 49 : 51. In 
Fig. 4, the proportion regulation under artificial partial extinction is shown. We start with 
the relaxed state of previous simulation (shown in Fig. 3(c)) with removing the 49 units 
which belongs to the cluster with negative u value (hatched in Fig. 4(a)). The remained 51 
units make a two clusters state with proportion p : 1 — p = 23 : 28 in the region B again. 
Three or more clusters state have not been observed. 

To clarify what selects the final state proportion, we perform numerical simulations of 
equations ([!]) with changing initial conditions. Parameters are the same as in Fig. 1 except 
K 2 = 0.5. We start from 500 initial conditions which are uniform random numbers between 
—0.1 and 0.1 with different seeds. Fig. 5 shows a distribution of finally selected proportion 
value with a peak at p = 0.5. These results show that the initial condition determines the 
proportion between the two clusters. 

Finally, we check a structural stability of these result by adding a small positive constant 
term e, i.e., 

itj = auj — bvj — if| +Ki(u — Uj) , 

Vj = cuj — dvj +K 2 (y — Vj) +e. 

Numerical simulation shows that the distribution of proportion also has a finite width. The 
most probable proportion, however, moves from 0.5 because the added term e breaks the 
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symmetry of ([[]) fL9| . For example, the most probable value is about 0.25, p m i n ~ 0.1, and 
Pmax ~ 0.4, respectively, when a = 0.6, b — 1, c = 2, d — 1, K\ — 0, K 2 — 3, and e = 0.2. 
From this result we conclude that the proportion regulation phenomena discussed in this 
letter is generic. 

The DIF which regulate the proportion of two types cells are widely studied about the 
differentiation process of Dd, e.g., cyclic adenosine 3'-monophosphate (cAMP), ammonia, 
and concentration of cation are known as candidates. If we specify the inhibitor and if 
we control the susceptibility, the proportion between two kind of cells is expected to be 
controlled. 

The authors wish to thank N.Nakagawa, Y.Kuramoto, and Y.Sawada for many fruitful 
discussions. The present work is supported in part by the Japanese Grant-in- Aid for Science 
Research Fund from the Ministry of Education, Science and Culture(No. 00060267). 
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FIGURES 

FIG. 1. (a) Temporal evolution of the distribution of units with respect to u. Gray scale 
represents the number of units by changing from to N between white and black. Randomly 
distributed units aggregate into the origin and one cluster state is formed, (b) Initial distribution 
of units. Each unit has a uniform random number between —0.1 and 0.1 for u and v, respectively, 
(c) Snapshot of one clusters state at T = 200. 

FIG. 2. Typical phase diagram of one and two clusters state. The upper line is given by 
(10) and the lower line denotes 1 — p c . The dotted line is K2 = K2 C . In the region A, there is 
no two clusters solution. In B and C, there is a linearly stable and unstable two clusters solution, 
respectively. 

FIG. 3. (a) Evolution of the distribution with respect to u. Units around the origin separate 
into two clusters state, (b) Initial distribution. Each unit has a uniform random number between 
—0.01 and 0.01 for u and v, respectively, (c) Snapshot of two clusters state at T = 500. 

FIG. 4. (a) Evolution of the distribution. After a sudden decrease of u of each unit to zero till 
T ~ 10, 51 units separate into two clusters again, (b) Initial distribution. 49 units in the hatched 
cluster are removed, (c) Snapshot at T = 500. Proportion p : 1 — p = 23 : 28 is selected. 

FIG. 5. Probability distribution of selected proportion of two cluster state from 500 randomly 
chosen initial conditions. 
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